Disordered asymmetric simple exclusion process: mean-field treatment 
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We provide two complementary approaches to the treatment of disorder in a fundamental nonequi- 
librium model, the asymmetric simple exclusion process. Firstly, a mean-field steady state mapping 
is generalized to the disordered case, where it provides a mapping of probability distributions and 
demonstrates how disorder results in a new flat regime in the steady state current-density plot for 
periodic boundary conditions. This effect was earlier observed by Tripathy and Barma |jj but we 
provide treatment for more general distributions of disorder, including both numerical results and 
analytic expressions for the width 2Ac of the flat section. We then apply an argument based on 
moving shock fronts |3 to show how this leads to an increase in the high current region of the phase 
diagram for open boundary conditions. Secondly, we show how equivalent results can be obtained 
easily by taking the continuum limit of the problem and then using a disordered version of the well- 
known Cole-Hopf mapping to linearize the equation. Within this approach we show that adding 
disorder induces a localization transformation (verified by numerical scaling), and Ac maps to an 
inverse localization length, helping to give a new physical interpretation to the problem. 

PACS numbers: 05.60.-k, 05.50.+q, 05.40.-a, 64.60.-i 



I. INTRODUCTION 



Many "real-life" nonequilibrium situations contain 
some kind of randomness / disorder and empirical observa- 
tions for example in traffic flow illustrate that such disor- 
der can lead to interesting new phenomena. In nonequi- 
librium statistical mechanics even one-dimensional (ID) 
models can exhibit phase transitions (see e.g., the re- 
view by Evans |3() and we are particularly interested 
in the effects of disorder on these transitions. Studies 
based on simple lattice based exclusion models incorpo- 
rate collective effects while offering possibilities for ana- 
lytic progress and easy computer simulation. 

In this paper we concentrate on the effect of quenched 
substitutional disorder on one such lattice model — the 
well-known asymmetric simple exclusion model (ASEP). 
The ASEP is one of the simplest nonequilibrium mod- 
els with a boundary-driven steady state phase transition 
and thus plays a paradigmatic role in nonequilibrium sta- 
tistical mechanics much as the Ising model does in the 
study of equilibrium systems. The present work is en- 
tirely within the framework of a mean-held approxima- 
tion and largely for the steady state but already shows 
many interesting effects on the phase transition such as 
an altered phase diagram and the presence of "Griffiths 
phases" 0. We hope subsequently to extend and com- 
pare this study with treatments allowing for fluctuation 
effects. Previous approaches to disorder in the ASEP, 
and related models, can be found in the work of Krug 
Kolwankar and others. Furthermore, field-theoretic 
approaches which retain the fluctuations can be applied 
to higher-dimensional generalizations of the continuum 
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FIG. f: DASEP with open boundary conditions. Filled circles 
denote particles, open circles are vacancies. The hard-core 
exclusion rule means the model incorporates collective effects. 
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version of the ASEP 

The paper is organized as follows. In Section [n] we 
define the model and summarize relevant results for the 
pure case. In Section IlIII we outline our two main meth- 
ods: a steady state mapping and a disordered gener- 
alization of the Cole-Hopf transformation. These two 
approaches are then developed further in Sections 
and IVII respectively allowing us to characterize (both 
quantitatively and qualitatively) the effects of disorder. 
In Section[V]we show how this affects the current-density 
diagram for periodic boundary conditions and the phase 
diagram for open boundary conditions and discuss finite 
size effects. Finally, in Section IV 111 we summarize our 
results and discuss areas for future work. 



II. THE DISORDERED ASYMMETRIC SIMPLE 
EXCLUSION PROCESS 

A. Definition of model 

The general form of the bond-disordered asymmetric 
simple exclusion process (DASEP) is summarized by the 
schematic of Fig. ^ A particle at site I hops to a va- 
cant nearest neighbor site on the right (left) with rate 
Pi (qi—i)- I n a discrete time version of the model (as 
implemented in simulations), these rates are replaced by 
probabilities per time step and a random sequential up- 



2 



i 




FIG. 2: Fundamental diagram (particle current J versus den- 
sity g) for pure ASEP with periodic boundary conditions in 
the thermodynamic limit. Note the particle-hole duality. 



date rule is applied. Here we consider for simplicity only 
the totally asymmetric case, qi — 0; qualitatively similar 
results are expected in the partially asymmetric case. 
Two obvious choices of boundary conditions are 

• Periodic boundary conditions: A particle from site 
L can hop into a vacancy at site 1 with rate pi,. 

• Open boundary conditions: Particles are injected 
at the left-hand end of the lattice and extracted 
at the right-hand end, forcing a particle current 
through the system. The usual convention is to 
insert particles onto site 1 with rate a if the site is 
empty and remove particles occupying site L with 
rate (3. It is this case which is illustrated in Fig. ^ 

Another possibility (used for example in the work 
by Popkov and Schiitz 0) is to have fixed "reservoir 
densities" g~ and g + respectively at left and right 
ends of the chain. By current conservation at the 
boundaries we see that the correspondence between 
reservoir densities and input and output rates a 
and P is given by 



a=p g (3=pl(1-Q + ) 



(1) 



where p$ is the hopping rate from the reservoir site 
to the first proper site 1 and pi, the rate from the 
last site L to the reservoir at the right. 

In the pure case fixing reservoir densities is exactly 
equivalent to fixing input and output rates, but in 
the case where the p/'s are disordered the two defi- 
nitions are different. We find both types of realiza- 
tion in representative problems e.g., traffic flow. 



B. Summary of results for pure case 
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FIG. 3: Phase diagram for pure ASEP with open boundary 
conditions. The second-order transition between high current 
(a/p > 1/2, j3/p > 1/2) and low current phases is represented 
by a solid line, while the first-order transition is shown dashed. 
For an infinite system the current in the high current phase 
is p/4; for a finite size system a slightly larger current can 
be sustained. The insets show sketches of the typical density 
profiles (site density gi against lattice site /) in each region. 



quantum-spin system |12| . Recall that the exact so- 
lution for the steady state can be summarized by a 
current-density plot (the "fundamental diagram" ) for pe- 
riodic boundary conditions and a phase diagram for open 
boundary conditions. For later reference these are shown 
in Figs. |3 and respectively. 

A mean-field approximation reproduces the exact 
phase diagram and also gives the essence of the dynam- 
ics through a treatment based on moving shock fronts. 
However, as might be expected from the low dimension- 
ality, it gives incorrect values for the two static exponents 
and the dynamic exponent. Similarly, we expect that a 
mean-field approach to the disordered problem will elu- 
cidate crucial features which are not controlled by fluc- 
tuations, such as modifications to the phase diagram and 
the dynamics but not the critical behavior. 

Recent work by Enaud and Derrida [13J has looked at 
the affect of disorder on the first-order phase transition, 
in this paper we concentrate mainly on the affect of disor- 
der on the second-order transition between low and high 
current phases. 



III. TWO PARALLEL APPROACHES 

In this section we outline our two basic approaches, 
explaining their use in the pure case and indicating how 
we extend them to treat disorder. 



The pure ASEP with site-independent hopping rates 
(i.e., pi = p for all I) can be treated by a variety of ap- 
proaches including recursive techniques Q, the ste ady 
state operator algebra formalism of Derrida et al. [T(| 
and its dynamic generalization [Tlj . or mapping to a 



A. Steady state mapping 

For general quenched substitutional disorder we can 
use our knowledge of the hopping rules to write an ex- 
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act expression for the average current across the bond 
between sites I and 1+1: 

Ji,i+i = (pi m(l - ni+%)) (2) 

where n\ = 1 or for a particle of vacancy at site I and 
the angular brackets denote an average over histories for a 
fixed pi. We introduce the average density for each site, 
(ni) = gi, and in the mean-field approximation ignore 
correlations between sites so 

(nm+i) (m)(ni + i} = QiQi+i- (3) 

When the system has reached a steady state the densities 
are constant in time and hence from the continuity equa- 
tion the current must be constant in space i.e., Ju+i = J 
for all I. So, for the mean-field steady state we have 

J = PlQlO- - Ql+i) (4) 
which then gives a mapping for gi+\ in terms of gr. 



For a particular realization of p;'s, if we know J and one 
of the densities (say g% ) we can use this mapping to ob- 
tain the density profile for the whole system. As we shall 
discuss in detail below, J is limited by the requirement 
that for all I we must have < Qi < 1. 

Note that we can also rearrange (@J to give a mapping 
for Qi in terms of f?2+i- in terms of the hole density, 
(Ji = 1 — gi, this "backward" mapping is 

<n = l - — — ■ (6) 

Picri+l 

which has exactly the same form as equation (JSJ due to 
the particle- hole duality of the system. 

The mappings for the pure case where pi — p indepen- 
dent of / have been given previously [l4| . There one finds 
that for low currents, J < p/A, the mapping has two fixed 
points. Mapping in the direction of increasing / the fixed 
point with higher g is stable and the lower one is unsta- 
ble. It is clear from JBJ that for mapping in the opposite 
direction the stability of the fixed points is reversed and 
the high a (low g) one is stable. The resolution of this ap- 
parent paradox is simply that the steady state selection 
in a given case is determined by the boundary conditions. 
Among the possible profiles are "kink" type solutions of 
the steady state pure mean-field profile map, having the 
form 

gi = \ + h taml 4> tanh(Z0 + 0) (7) 

where tanh0 = — AJ/p) and 9 are determined by 
boundary conditions. At J — p/A the two fixed points 
combine in a half-stable fixed point at g — 1/2. In the 
high current regime, J > p/A, the mapping has no fixed 
points and density profiles have the form 

gt = i - I tan 4>' tan(Z<£' + 6') (8) 
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FIG. 4: Comparison of mean-field mapping and simulation 
data for densities in a system of 100 lattice points with open 
boundary conditions (pi's drawn from uniform distribution 
0.45-0.55, simulation with a = = 0.2). 



with tan <j>' = \J (AJ/p — 1). Since we must have < gi < 
1 for all / then JBJ clearly applies with non-zero <fr' (i.e., 
J > p/A) only if I is confined within the boundaries of 
a finite system. These mean-field profiles agree qualita- 
tively with exact solutions ^3] although the mean-field 
versions over exaggerate the sharpness of the shock front, 
which in practice is broadened by fluctuations. 

In the disordered case it is straightforward to iterate 
by computer. For specific realizations of disorder (i.e., 
particular choices of {pi}) we have compared densities 
from this mean-field mapping with profiles obtained by 
Monte Carlo simulation. As shown in the low current 
example of Fig. 0] there is a reasonable qualitative fit 
though again the mean-field shock front is sharper than 
the simulation result. Close to the pure critical current 
one sees that in some regions (corresponding to groups 
of "weak" bonds with low pi) the profile has roughly the 
high current form while in other regions ("strong" bonds) 
it takes the low current form — this is the analogue of Grif- 
fith's phases in a magnet fj| . We performed comparisons 
for the open boundary case (see |lj for similar discussion 
with periodic boundary conditions) where the simulation 
is controlled by parameters a and /?. A computer fitting 
procedure was used to match a and /3 to the value of J 
and one of the gis needed to implement the mean-field 
mapping (for the pure case this can be done exactly). 

Such examples for specific realizations of disorder sup- 
port the qualitative validity of the mean-field approach. 
In Section ll VI we develop a more general analysis to pre- 
dict the typical effects of a given distribution of disorder 
(Gaussian, uniform etc.). By considering the mapping of 
density distributions we are able to show that a charac- 
teristic effect of disorder is a shift in the average density. 
The physical meaning is explored further in Section Ivl 
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B. Continuum limit and Cole— Hopf transformation 

Here we consider the continuum limit of the DASEP 
and introduce a disordered generalization of the well- 
known Cole-Hopf transformation 0, 0] . 

For "smooth enough" disorder (and working once again 
in the mean-field approximation) we can take the contin- 
uum limit of equation J2J), to arrive at 



J = ^[e{i-Q)-~ 



(9) 



where J and g are in general functions of continuous po- 
sition x and time t, and we set the lattice spacing equal 
to 1 for convenience. Substituting this into the continu- 
ity equation yields a (noiseless) disordered Burger's type 
equation 



dg 
dt 



d_ 

dx 



p(x) 



,, s 1 dg 
e{l - Q) -2d-x 



(10) 



The next step is to transform to a height variable h, such 
that dh/dx = g — 1/2, giving 



d 2 h 
dxdt 



d_ 

dx 



p(x) 



dhy l d 2 h 
dx J 2 dx 2 



(11) 



This can be trivially integrated with respect to x to give 
a noiseless disordered version of the gro wth model stud- 
ied by Kardar, Parisi and Zhang [l7|. We then put 
h(x, t) = A ln[u(a;, t)]+f(t) and choose the arbitrary func- 
tion f(t) and the constant A so as to remove all non-linear 
terms. This disordered generalization of the Cole-Hopf 
transformation finally gives us 



du 

at 



D(x) 



d 2 u 
dx 2 



D{x)u 



(12) 



where D(x) = p(x)/2. Equation (|12|l is a linear equa- 
tion (therefore much easier to treat numerically and an- 
alytically) which still preserves the full dynamics of the 
system. From its solution g is given by the inverse Cole- 
Hopf transformation 



Q - 



1 1 d In u 

2 = 2 dx 



(13) 



For a steady state solution for g then u must be separable 
with time dependence e~ ut and lu = 2J. 

In the pure case the second term on the right-hand side 
of 112|) can be absorbed into the definition of the trans- 
formation to leave us with just the diffusion equation, 



du j^p 2u 
dt dx 2 



(14) 



This diffusion equation can be trivially solved and map- 
ping back through the Cole-Hopf transformation then 
gives the well known pure continuum steady state solu- 
tions which are the continuum versions of J7J and © . 



For the disordered case, the situation is more compli- 
cated and in particular we have to include the —D(x)u 
term on the right-hand side of l|12(l . It is not immediately 
obvious how to treat this equation for general D(x), al- 
though it might be possible to solve it for specific D(x) 
or to do some kind of WKB type approximation. In Sec- 
tion ^T] we shall show how much useful information can 
be obtained via a more powerful numerical scaling ap- 
proach which reveals that disorder induces a localization 
transition. This provides a complementary approach to 
the steady state mapping outlined above and enables us 
to interpret the effect of disorder as a localization tran- 
sition in the transformed system. 



IV. STEADY STATE MEAN-FIELD MAPPING 
FOR DISORDERED CASE 

A. Mapping of distributions 

Here we return to look in detail at the steady state 
mapping with a known distribution of disorder. It is 
convenient to take the disordered variable as 



7i = J/ Pi- 



(15) 



If the position-independent distribution f(p/i) is known, 
then one can use the mapping to relate the probability 
distribution w of Qi+\ to the distribution of gf. 



Wl+i{Ql + i) 



Wl 



7 



11 



i - ei+ij (i - Qi+i) 



rf(l) dry. 



(16) 

The subscripts on the w's in Qlfijl indicate that we expect 
the distribution to change as we map through the system. 
For example, if we start from a known g\ (i.e., w% is a 
delta function) then the width of the distribution will 
obviously increase as we look at w 2 , etc. 

After iterating the mapping for many steps the density 
distribution will eventually converge on some fixed point 
shape, the position of which will depend on the direc- 
tion in which we map (just as in the pure case discussion 
in IIII A|) . This stationary probability density is the dis- 
tribution of Qi's which we would expect to see in the pe- 
riodic boundary case. However, as we shall demonstrate 
in IV Bl by considering moving shock type solutions we 
can also gain some information about the expected open 
boundary profiles and phase diagram. 

Numerically therefore, we look for the stationary prob- 
ability density of gi's for different J (averaged over many 
realizations of disorder). In practice this involves using 
the mapping of equation for large system sizes (say 
10,000 lattice points), repeating for different go and dif- 
ferent realizations of disorder then creating a histogram 
of the values of g. 

We find that a convenient order parameter to charac- 
terize the distributions is the asymmetry about g = 1/2 
given by A = (g — 1/2) (where numerically we take the 
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FIG. 5: Numerical mapping data for A = {g — 1/2} versus 
J in pure case and specimen disordered cases. Pure case has 
fj = 2.0 {rj denned as 1/p) corresponding to a pure critical cur- 
rent Jq of 0.125, disordered cases are Gaussian distributions 
with the same mean and standard deviations a v = 0.2, 0.5. 
Lines are provided as an aid to the eye. For arbitrarily weak 
disorder, A is non-zero for all J. 



average over all g in the physically accessible region be- 
tween and 1). In the pure case we find that A is zero in 
the high current regime and non-zero in the low current 
regime (see Fig. |3J). This is just what we would expect 
from using the pure mapping (see again ITTl Afl in the for- 
ward direction: in the low current phase the densities 
in the bulk of the system will approach the upper fixed 
point, whereas for a long system in the high current phase 
the densities in the bulk are all very close to 1/2. How- 
ever, in the disordered case we find that A is non-zero 
for all J as shown for typical examples in Fig. [5] 

We now consider an analytical argument to reproduce 
this data basing our method on the work of Hirota 0] 
for different random mapping processes. In our case, we 
impose the fixed point condition by setting wi — = 
w so that llt)|) gives (dropping redundant subscripts) 



l- e )(l- e ) 



:f(l)dr (17) 



This integral equation is difficult to solve analytically 
but we can find approximate solutions by considering the 
dominant terms in different regimes. In IIV Bl ITvTjl we 
pursue this approach above, below and close to the pure 
critical point. Then in IIV El we consider the specific sol- 
uble example of a Lorentzian distribution of disorder. 



an expansion of the integrand in powers of a 2 and obtain 
a functional differential equation 



w(g) = w 



2 dg 2 



l-Q 



(18) 



where terms involving higher moments have been ne- 
glected. We now wish to solve (11811 for normalized non- 
negative w(g). Note that in order to be able to carry out 
the integrals analytically we here allow g to take any real 
value whereas in the physical problem < g < 1. As we 
shall discuss later this is not expected to introduce much 
of an error providing J is not too high compared with 
the pure critical current. For a 2 "small" (in a sense to 
be clarified) we assume a solution of the form 



w(g) = w (g) + 

The pure solution Wo(g) is easily shown to be 

1 



w (o) = A- 



Q ~ Q + 7 



(19) 



(20) 



with A a normalization constant given by y/'y — 1/4/tt. 
This solution is only valid for 7 > 1/4 which corresponds 
to being in the high current phase of the corresponding 
pure model. Note that we can also get this result from 
the known high current mean-field pure result © using 
the obvious relationship w(g) ~ l/|jf |. By considering 
this form for w (g), we see that to satisfy l(T%|) Wi(g) 
must have a factor (g 2 — g + 7) 3 in the denominator. 
Straightforward calculation gives 



wi(g) 



2A 



Q - 1Q 



1 {q 2 - Q + lY 



(21) 



From this expression for the stationary probability dis- 
tribution we can calculate (via contour integration) the 
average value of g — 1/2: 



(q- \) w(g) dg 



4(7 -i)' 



(22) 
(23) 



This is the chief analytical result of this subsection; we 
now compare it with numerics and discuss its validity. 

For comparison with data it is more helpful to 
write l|23(l in terms of the current J and the mean fj and 
variance a 2 of the inverse hopping probability r\i = 1/pi, 



B. Calculation of A above pure critical point 

Above the pure critical point the integral in (|17fl is 
dominated by the peak in 7(7). For disorder sharply 
peaked about 7 with small variance a 2 we can perform 



2^2 



A 



Pa 



(24) 



In Fig. H3 we compare the prediction of this analytical 
result with the data from our numerical mapping with 
a Gaussian distribution of disorder. We see that l|24|) 
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FIG. 6: Comparison of numerical and analytical results for 
A, above pure critical current {J% = 0.125). Numerical data 
for Gaussian disorder with t) — 2.0 and a v = 0.2. For weaker 
disorder the analytic expression is valid for currents closer to 
the pure critical point — see discussion in text. 



reproduces well the general trend of the data but there 
are a couple of obvious problems. For high values of 
J, the analytical expression slightly under-estimates the 
numerical result — this is due to the fact that analyti- 
cally — oo < g < oo whereas in the numerics (as in the 
physical problem) we have averaged over < g < 1. 
As expected, this discrepancy becomes larger for high J 
since then the density profile is steeper and is concen- 
trated less about g ~ 1/2 so the tails of the distribution 
are more important. One can numerically integrate (e.g., 
using maple or similar) the expression in (12211 between 
and 1 (and adjust the normalization correspondingly) 
and obtain values in better agreement with the data. 



A more serious problem is the fact that our analytical 
expression for A tends to infinity as J tends to the pure 
critical current Jq = 1/ (Aff) , whereas numerically we see 
no divergence in A. In fact it is easy to see why the 
analytical method fails close to the critical point. As 
J -> j£, the pure distribution wq{q) becomes more and 
more sharply peaked about g = 1/2, however W\(g) will 
be even more sharply peaked (due to the cubed term in 
the denominator) and eventually the magnitude of this 
first order correction becomes large compared with wo(g) 
for some values of g. This leads to unphysical results such 
as negative probability densities. A rough calculation 
shows that in order for the perturbation expansion of l|18|) 
and l(TT)|) to be valid we require a 2 <C (7 — 1/4) and 
this condition becomes impossible to satisfy as 7 — > 1/4. 
Alternative approaches to obtain an expression valid near 
the critical point will be considered in II V Dl below. 



C. A below pure critical point 

Far below the critical point the integral in l|17l) is 
dominated by the contribution from the sharp peak of 
w (jzr^)- Performing a saddle-point expansion about this 
peak (and for convenience assuming a Gaussian distribu- 
tion for the disorder) we obtain 



g[w(g)} 2 exp(- 



>{Q) = 



2^H 



a[w"(g)w{g) - [w'(g)} 2 



1/2- 



(25) 



Now in the saddle-point expansion g is defined as the 
maximum of w(g) so for consistency it must correspond 
to the maximum of the Gaussian in l|25|) . i.e., we have 



(1 - g)g~ = 7 



(26) 



which is just the condition for the fixed-points of the pure 
low current mapping. We recall that the upper of these 
two fixed points is stable as we map forward through the 



system while the lower one is unstable. So as cr^ 



0, 



w(g) of equation (|25|l tends to a delta function about the 
upper fixed point. The addition of disorder, broadens 
this pure distribution to a Gaussian with the same mean 
g and standard deviation cr^/g, i.e., 



;(<?) = (2^)- 1 l 2 {~g/a 1 ) exp 



(q-q) 2 



(27) 



2{a 1 /gf 

Calculating A in this approximation is trivial and gives 



A = 



7- 



(28) 



So we conclude that, for currents well below the pure 
critical current, disorder doesn't change the pure result 
for A. This agrees with the numerical results for J <C Jq 
shown in Fig. 

Again our analytic prediction fails close to the critical 
point; this is due to the breakdown of the assumption 
that the major contribution to the integral equation 117|) 
is due to the sharp peak in w(g). From the above analysis 
we can see that the standard deviation in 7 of the peak 
in ui(t^-) is about er 7 (l — g)/g so if this peak is to be 
sharper than the one in /(t) we require: 



<7 7 (1 — g)/g <C cr. 



g > 



(29) 
(30) 



This condition breaks down as we approach the pure crit- 
ical point and a more sophisticated analysis becomes nec- 
essary (see Subsection IIV Dt 

Note that in order to obtain an analytic approximation 
for A in the low current regime we needed to assume a 
particular distribution for the disorder. Repeating the 
procedure for different distributions (e.g., uniform, bi- 
nary) we find that w(g) has a different form in each case 
but to a first approximation A takes the pure value in 
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FIG. 7: Comparison of numerical and analytical results for 
A, below pure critical current (J c = 0.125). Numerical data 
for Gaussian disorder with 77 = 2.0 and a v = 0.2. Note the 
discrepancy close to J c . 



of disorder where the mapping equations turn out to be 
exactly soluble. This particular case helps build up a 
general picture of what happens near the critical point. 



E. An exactly soluble case 

Hirota and Ishii have treated exactly the case 
where the disorder variable in their mapping has a 
Lorentzian distribution. They show that the stationary 
probability density is also Lorentzian and calculate its 
width and mean. Our mapping equation has a different 
form but is also amenable to analytical treatment for the 
case of Lorentzian disorder. This is therefore a useful 
test case where we can compare numerics and analytics 
for the full range of J. 

If 7 is drawn from the Lorentzian distribution 



/(7) - - 



r 



7r (7 — 7) 2 + r 2 ' 



(33) 



all cases. This is in contrast to the high current regime 
treated in llV Bl where the exact form of disorder is irrele- 
vant to first order (the relevant parameter is the standard 
deviation of 1/p) and w(g) takes the universal form given 
by (|19l) - (|21|l . The numerics confirms these arguments. 



D. Calculation of A around pure critical point 

Ideally we would like an expression for A at and very 
close to the pure critical current since this is the most 
interesting regime physically (a large system cannot sus- 
tain currents much above the critical point as will be 
discussed in detail in Section below) . Unfortunately 
in this intermediate regime it is not easy to see how to 
treat l|17|) as both w(j^) and f(j) are sharply peaked. 

One approach is to assume that the product 
w( j^—)f{ r y) is sharply peaked in 7 and perform a saddle- 
point expansion on this. For convenience, we define 
y{g) = — In and consider a Gaussian distribution 

of disorder. Then the saddle-point x is defined by 



= 



[x(l - q) - 7] _ _ 1 



(l-g)--+y'(x) (31) 



where the prime denotes differentiation with respect to 
the argument of the function. Evaluating l|17|) about this 
saddle-point yields 



y(e) = 



[x(l - g) - 7f 
2a 2 



-hxx + y(x). (32) 



In principle these equations are sufficient to determine x 
and y(g) but the implicit definition of x makes further an- 
alytic progress extremely difficult. Instead we consider, 
in the next section, a specific (Lorentzian) distribution 



and we assume a Lorentzian distribution for gi (mean 
g, width t) then we can integrate exactly the integral 
mapping equation Ijl7|l to find that gi+i also obeys a 
Lorentzian distribution with mean and width given by: 



1 



1Q- 



Tt 



t' = 



g 2 +t 2 
Tg + ^ft 
p 2 +t 2 ' 



(34) 
(35) 



Following ^3 we can characterize each Lorentzian dis- 
tribution by a complex number whose real part represents 
the mean and whose imaginary part represents the the 
width. Then the mapping relationship is 



R' = 1 - Gfl/R 



(36) 



where R — g + it, G — 7 + iT and the dagger denotes 
complex conjugation. Now the fixed point of the mapping 
is given by the Lorentzian distribution characterized by 
R* = g* +it*, with 



Q = 
t 



t* + r 

2t* 



= x 1 ' 2 



where x is the positive root of 

x 2 + a - 7 ) x - 



0. 



(37) 
(38) 

(39) 



From this stationary distribution one can calculate A; 
averaging over g's from — 00 to +00 gives 



A 



Wx- 1 / 2 . 



(40) 



As discussed above (see text following equation l|24|) ). 
this extension of the range of g is not expected to make 
much difference unless the current J is large. In fact for 
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FIG. 8: Comparison of numerical data and analytical pre- 
dictions for Lorentzian distribution of disorder with fj = 2.0 
and width=0.2. Relatively large width means affect of r\ < 
cutoff is significant. 



this Lorentzian case we can actually do the integral over 
g between and 1 analytically, yielding the proper result 



1/2 / ^1/2 _ r) 2 + 4j ,2 



2tt 



In 



(a; 1 / 2 - 

r 



r) 2 



■ 4a; 2 



2ttx 1 /2 



arctan 



4a .3/2 



r 2 + 4a; 2 



(41) 



In Fig. we compare the results of these expressions with 
numerical data for a Lorentzian distribution of disorder 
and find as expected that l|41|) gives a noticeably better 
fit than l|4U|) for J above the critical point. The fit is still 
not exact due to the fact that in the numerics we have 
imposed the physical restriction that 7 must be positive 
but this cutoff is not incorporated in the analytics. This 
problem is more pronounced than in the Gaussian case 
because of the relatively high weight in the tails of the 
Lorentzian distribution. 



V. EFFECT OF DISORDER ON 
FUNDAMENTAL AND PHASE DIAGRAMS 

We have seen that the shift A provides a useful char- 
acterization of the effect of disorder; we now turn our 
attention to what this shift means physically in terms 
of the density profiles, fundamental diagram (for peri- 
odic boundary conditions) and phase diagram (for open 
boundary conditions). In Section llVI we presented both 
numerical and analytical approaches to calculate A for 
fixed J. However in simulations J is not held constant 
and we must consider how it changes when we add disor- 
der. We concentrate initially on the thermodynamic limit 
where J is easier to predict and then in IV CI consider the 
importance of finite size effects in small systems. 



Recall that the treatment of Section HVl was based on 
using equation JSJ) to map forwards in g. If instead we 
had used © to map backwards in a then by exactly the 
same argument we would have concluded that the steady 
state distribution for the disordered case has (<r) > 1/2 
i.e., (g) < 1/2 for all J. This apparent paradox is exactly 
analogous to the situation in the pure low current phase 
where the stable fixed point value of g depends on which 
direction we map (see discussion in llll Afl . The resolution 
in both cases is that which fixed point is seen in the bulk 
depends on the boundary conditions. We now consider 
these in more detail. 



A. Periodic boundary conditions 

Perhaps the most obvious change expected when we 
add disorder to a large system is a decrease in the max- 
imum sustainable current. For an infinitely large pure 
system with periodic boundary conditions the maximum 
current J^ axp is just the critical current J®, = p/4. 
Similarly for a disordered system, the maximum pos- 
sible current is limited by stretches of "weak" bonds 
(i.e., low p) so in the thermodynamic limit we expect 



J maXiP = Jc = Pmin/4 (where p m in is the smallest value 
of p permitted by our distribution of disorder). Using 
the methods of Section ITVl we can obtain A correspond- 
ing to all possible currents up to Jc- In contrast to the 
pure case we now have a non-zero value of A( Jc) which 
we shall denote for convenience by Ac- We now address 
how this is reflected in the fundamental diagram. 

Let us first consider currents below the maximum. 
Clearly, one possibility is for the distribution of g for 
all sites in the lattice to be given by the stationary dis- 
tribution w(g) determined from the condition l|17f) . The 
average density is then obviously given by 1/2 + A. How- 
ever, the argument above illustrates that it's also possi- 
ble for the distribution to be at the unstable fixed point 
of the forward density mapping giving (a — 1/2) = A 
and hence an average density of 1/2 — A. So, in the 
low current phase with periodic boundary conditions the 
possible disordered profiles are roughly the same as in 
the pure case (i.e., either at the upper fixed point or the 
lower fixed point) although of course with added noise. 
The exact position of the fixed points (characterized by 
A) deviates slightly from the pure case especially close 
to the critical point. 

However, in the maximum current phase it is possible 
for the distribution to start near the lower unstable fixed 
point and map forward to the upper stable fixed point via 
a noisy shock front (whose position may alter). The pe- 
riodic boundary conditions are maintained by stretches 
of decreasing g corresponding to weak bonds in the pure 
high current phase. A crude way to consider this is to 
look at the profile as a superposition of a high current 
profile with a small shock front type low current profile 
as shown schematically in Fig. [5] This "density segre- 
gation" into sections (not necessarily of equal lengths) 
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FIG. 9: Possible density profile in maximum current phase for 
system with disorder. Periodic boundary conditions can be 
imposed for example between points BC, AF or DE leading 
to a range of possible densities from 1/2 — Ac to 1/2 + Ac- 




FIG. 10: Schematic fundamental diagram for DASEP in ther- 
modynamic limit. 



with density 1/2 — Ac and 1/2 + Ac) was previously 
explained for the binary disorder case by Tripathy and 
Barma [[J . Enforcement of the periodic boundary condi- 
tions (see again Fig. then leads to macroscopic average 
densities in this maximal current phase anywhere from 
1/2- Ac to 1/2 + Ac). 

So the end result is a fundamental diagram which looks 
like the pure one for low currents but has a new flat 
regime of width 2Ac at the maximum current as shown 
in Fig. ^3 This flattening effect was observed in Q] and 
the width of the flat section calculated for the particu- 
larly simple case of binary disorder. Our method enables 
us to treat more general distributions of disorder — from 
A(J) and the maximum current then we can construct 
the complete mean-field fundamental diagram. 

Figure 1111 shows the fundamental diagram obtained 
from Monte Carlo simulations for a system of size 5000 
with a particular realization drawn from a uniform distri- 
bution of disorder (results for the pure case corresponding 
to the mean of 1/p are also shown) . Comparison with the 
calculated A( J) for the same width of disorder (Fig. [SJ) 
shows that the mean- field prediction is roughly correct, 
though the maximum current in the simulation is larger 
than the mean- field thermodynamic prediction of 0.1. 
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FIG. 11: Monte Carlo simulation data for DASEP fundamen- 
tal diagram for system size 5000. The uniform disorder case 
has fj = 2.0 and a v = 0.5; flattening with respect to the pure 
case is clearly seen. 



B. Open boundary conditions 

Popkov and Schiitz have shown how to predict the 
phase diagram for open boundary conditions from the 
fundamental diagram for periodic boundary conditions. 
Their argument considers the motion of shock fronts 
thorough the bulk to motivate an extremal current prin- 
ciple: 

J = max J(g) for q~ > g + (42) 
ee[e+.e _ ] 

J = min J(g) for g~ < g + (43) 
ee[e-,e+] 

where g~ and g + are reservoir densities fsee III AT . 

Applying this Popkov-Schiitz argument to the funda- 
mental diagram of IV Al leads to a growth in the size of 
the high current phase (as compared to the pure case) 
resulting from the flat section on the fundamental dia- 
gram. This growth is by an amount Ac in both the g~ 
and the g + direction. The resulting phase diagram in the 
g~—g + plane and comparison with the pure case is shown 
schematically in Fig. 1121 So our numerical and analytical 
calculations of A (J) in the previous section allow us to 
determine quantitatively the disordered phase diagram; 
compare the simulation data of Fig. ED 



C. Finite Size Effects 

In this subsection we outline briefly the modifications 
to the above picture for finite size systems. The discus- 
sion is inevitably fairly qualitative and it is worth noting 
that even in the pure case, a mean-field treatment does 
not correctly capture all finite size effects. 

In the pure case there are two main finite size effects. 
Firstly in the open boundary case the system can sustain 
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FIG. 12: Phase diagram for DASEP in the g~-g + plane (black 
lines) with pure case for comparison (grey lines). 




FIG. 13: (Color online). Monte Carlo simulation results for 
DASEP phase diagram — surface and contour plots showing 
current as a function of g~ , g + . Uniform distribution with 
f\ = 2.0 and a v = 0.5, system size 5000. Note increase in size 
of flat maximum current area compared to the pure case. 



field theory which predicts J c = p/4 for all system sizes. 

In the disordered case Jc is limited by stretches of 
weak bonds and the probability of finding within the sys- 
tem a long stretch of consecutive weak bonds increases 
with system size. Hence in mean-field theory we expect 
Jc to be greater than p m in/4 but less than the corre- 
sponding pure result. It is relatively straightforward to 
calculate the expected mean-field Jc for simple cases 
such as a binary distribution (see e.g., Q). However, 
just as in the pure case, we expect the true value of Jc 
for small systems to be larger than this mean-field es- 
timate. In addition, a novel feature of the disorder is 
that it allows density profiles in the periodic boundary 
case such as that shown in Fig. where the current flow 
is larger than Jc, i.e., for finite disordered systems we 
can have J maXjP > Jc- It is even possible to conceive 
of situations in which J ma x,p is increased above J c by 
adding disorder, meaning that A above the pure critical 
point (as calculated in IIVB|I can be a physically rele- 
vant quantity. Examination of the possible high current 
profiles reveals that for small systems the current varies 
with macroscopic density (with the maximum J ma x,p at 
g = 0.5) so a flat section is not seen in the fundamen- 
tal diagram. Similarly in the disordered finite size open 
boundary case, J max ,o is increased beyond Jc and there 
are alterations in the phase diagram corresponding to the 
altered fundamental diagram. 

So finite size effects have a significant complicating in- 
fluence on both the fundamental diagram and the phase 
diagram. At present we are not able to quantify these en- 
tirely even within mean-field but progress can be made by 
combining numerical work (e.g., self-consistently looking 
for the maximum current the density mapping can sus- 
tain) with the analysis of previous sections. 



VI. RESULTS FROM DISORDERED 
COLE-HOPF TRANSFORMATION 

Recall from Subsection IIII Bl that we are able to treat 
the continuum limit of the DASEP via a disordered gen- 
eralization of the Cole-Hopf transformation (|13|l to ob- 
tain the linear equation (|12|l . Here we develop this ap- 
proach further and demonstrate connections to the re- 
sults from the steady state mapping. 



A. Scaling and localization in the steady state 



a current J max , which is slightly larger than J c . One can 
obtain a mean-field prediction for </ max . by looking for 
the largest value of J for which all gi of © are in the 
physically applicable regime between and 1. Secondly, 
in a small system J c is increased slightly from p/A — 
this effect is due correlations between particle densities at 
adjacent sites and is therefore not reproduced by mean- 



Let us concentrate initially on the steady state solu- 
tion for g in order to make comparisons with the discrete 
mean-field mapping approach. As discussed in IIII Bl a 
steady state solution for g corresponds to a separable 
solution for u i.e., u(x) — T(t)X(x). The a>dependent 
factor X must then satisfy 

-[u-D(x)]X=D{x)^ (44) 
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with to = 2 J. 

Our approach is to rediscretize this, 

(3D n — u>)X n = D n (X n+ i + X n _i), 



(45) 



and then employ numerical scaling based on a method 
developed by Pimentel and Stinchcombe [2(| to treat the 
equation of motion of a ID Mattis-transformed Edwards- 
Anderson Heisenberg spin glass. We write l|45[l as 



i X„ 



V nn +iX. 



n+l 



(46) 



with £„ = D/D n , E n — 3D, V n , n +\ = D where D is the 
characteristic strength of the disorder variable. Equa- 
tion i|46|) is of just the form considered in [2(| and can be 
exactly scaled by b = 2 decimation to give 



(E' n - ( n w)X r , 



= V 

n.n- 



2 A„_ 2 + V nn+2 X n 



+ 2 



with 



T/ / _ Ki,n+l^n+l,ri+2 



E 



n+l 



E„ — E„ 



r2 

n.n — 1 



V 2 



V 2 ^ 



(47) 

(48) 
(49) 



So the E n and V„ !Tl ±i, which are n-independent at the 
outset, pick up correlated randomness under scaling. 

It is easy to iterate these equations numerically and 
check how the V's evolve. For the pure case (i.e., all 
Cn = 1) we find an allowed energy band for D < to < 5D. 
Within this band the "site potential" E and the "cou- 
pling" V evolve chaotically, corresponding to extended 
states. Outside the band there are no allowed states and 
V decreases rapidly and monotonically to zero while E 
tends to a constant value. Just as in [2(| this can be ex- 
plained by writing V and E explicitly in terms of a single 
parameter 9 which is related to the wave vector of exci- 
tations within the band. As expected the lower edge of 
the allowed band lu = D, corresponds to the pure critical 
current Jq = p/A. The upper band edge has no physical 
significance in our problem since this switch between con- 
tinuum and discrete representations is valid only for long 
wavelengths corresponding to being close to u> = p/4. 

Adding weak randomness, we find that for all values 
of u, V evolves either chaotically or cyclically to zero in- 
side a well defined exponentially decreasing envelope i.e., 
V{1) ~ f(l)e~ L /£ where I is the distance between sites and 
£ is a localization length. In other words, any amount of 
disorder induces localization for all values of frequency. 
This is analogous to the fact that all states are local- 
ized in one-dimensional disordered quantum problems. 
We developed a computer algorithm to calculate £ for a 
given distribution of disorder (averaging over many real- 
izations) as a function of to. Typical results for both pure 
and disordered cases are shown in Fig. ^31 Note that, in 
contrast to the spin-glass case studied in |2jj, the local- 
ization length is not infinite at the critical point. This 
is essentially due to the —D(x)u term in l|12f> and means 
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FIG. 14: Localization transition in Cole-Hopf mapped 
DASEP. The data points represent the localization length £ 
calculated by our numerical scaling method; lines are provided 
as an aid to the eye. In the pure case (with fj = 2.0) the lo- 
calization length is effectively infinite in the "allowed band" 
0.25 < u> < 1.25: disorder is seen to induce localization. 



that we cannot easily follow |20j| in defining a dynamic 
exponent via a relationship like u) tx (l/£) z . 

To determine the signature of localization in the orig- 
inal DASEP problem we consider a localized form of u, 



(50) 



Using i|13fl to invert the Cole-Hopf transformation gives 



e 2 2£ 



(51) 



i.e., a shift in the profile. So calculating the average 
localization length is just like calculating the average shift 
from q = 1/2 in the original problem. In other words 
the average localization length is just half the inverse 
of the quantity A defined in Section IIVI Of course, in 
general u will be some boundary-dependent combination 
of e +x ^ and e~ x ^ , corresponding to the fact that, in the 
high current phase, the observed shift in density can be 
anywhere between +A and —A i.e., the flat section on the 
fundamental diagram. This relief of localization in the 
inverse mapping is possibly connected with the work of 
Kopidakis and Aubry [2l| on the relief of localization by 
non-linearity in low-dimensional deterministic systems. 

In Fig.^]we explicitly compare the localization length 
obtained by this scaling method and (2A) _1 from the 
discrete mapping method. We find an excellent agree- 
ment for J near to the critical point where the discrete- 
continuum-discrete approximations are valid, but the 
comparison breaks down in the high J region which 
in any case is unphysical. Furthermore, by defining 
Y n+ i = X n+ i/X n we can cast (|45|l into exactly the form 
studied by Hirota 0] 



Y, 



n+l 



l/Y n 



(52) 
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FIG. 15: Comparison of localization length £ obtained by 
numerical scaling with (2A)" 1 from the mapping approach. 
Disorder with f\ = 2.0, a n — 0.2. Analytical predictions for £ 
and (2A) _1 in the high current regime are also shown. 



with 



(3 - ( n u). 



(53) 



We can then obtain the stationary probability distribu- 
tion w(Y) in analogy with the calculation of Section llVl 
and, following 18], define the localization length £ by 



1 1 f°° 

7=0 w(Y)lnY 2 dY 

S ^ J — CO 



leading eventually to 



[4 - (3 - 2^) 2 ] 



2u 2 <j 2 



(54) 



(55) 



The resulting prediction for the localization length of 
the high current phase agrees closely with that deter- 
mined by our numerical scaling method (see Fig. I15|) 
except near the critical point where the Hirota method 
breaks down (just as in the pre Cole-Hopf case). One 
subtlety involves the meaning of "averaging over realiza- 
tions" in finding the average localization length. In the 
original numerical scaling method of pp| | the averaging 
is over V ~ e^ 1 ^ 1 (where i labels the specific realization 
of disorder) ; in contrast the Hirota method takes the av- 
erage of 1/&. It is clear from (|51|l that averaging over 
1 /£j provides the definition of £ most directly comparable 
with A so this is the procedure adopted in the computer 
programs used to generate the data in FigureslT^andlTol 
Averaging over V instead gave less good agreement with 
the analytical expression especially in the center of the 
band where the localization length is less clearly defined. 

Finally we note that since the variable Y can physically 
take all values between — oo and +oo we expect better 
agreement in the numerical and analytical results for £ 
than the numerical and analytical results for A where the 



analytics for general distributions of disorder was unable 
to take account of the physical restriction < g < 1 
resulting in discrepancies for high values of J (see dis- 
cussion of IIVB|I . In addition, it is easier to get high 
quality data from the linear scaling computer algorithm 
(where we can measure accurately localization lengths up 
to ^10,000) than from the non-linear mapping. 



B. Effect of disorder on dynamics 

One advantage of the Cole-Hopf formalism is that it 
provides an easy route to discuss the influence of disorder 
on the dynamics of the system. Here we mention briefly 
some inferred effects. A general solution for u is made 
up of a superposition of separable solutions; in the pure 
case this leads to a general solution for g of the form, 



g(x) 



1 _ 1 Y. k A k ike ikx - Wkt 

2 ~~ 2 Y, k A k e lkx -^kt 



(56) 



where k can be positive or negative and is real in the 
high current phase and imaginary in the low current 
phase. The coefficients A k must be chosen so that g 
is real. From (|12|l the pure "dispersion relation" is 
uj k = D(l + k 2 ) but the /c-independent term will can- 
cel out in the numerator and denominator of (|56|l . The 
solutions hence have wavelike form in the high current 
phase and multiple soliton form in the low current phase. 
As time increases, transients die away leaving the steady 
state corresponding to the smallest value of u>. 

Disorder induced localization in the high current 
phase, is crudely like adding a small imaginary part 
in (k ~ l/£) to the real k (in l|56|l this gives steady 
state solutions for g like that in Fig. |!5J). The imagi- 
nary part of the resulting complex dispersion relation 



D(l + k 2 



2ikn) would be expected to lead 



to oscillations while the small decrease in the real part 
slows down the dynamics. Indeed, it is intuitively obvi- 
ous that adding disorder should slow down the approach 
to the steady state, since in a ID system the overall rate 
of hopping of the particles (and hence the speed at which 
the steady state is reached) will be limited by the bond 
with the smallest pi . This slow-down is confirmed by the 
Monte Carlo simulations of Fig. Within this frame- 
work a more quantitative analysis should be possible but 
would be complicated by boundary conditions and finite 
size effects (compare Section |Vj. 



VII. DISCUSSION AND OUTLOOK 

In this paper we have shown within a mean-field frame- 
work that one effect of quenched bond disorder on the 
ASEP is a flattening of the top of the steady state 
current-density relation and a corresponding increase in 
the high current region of the phase diagram for open 
boundary conditions. We have presented various numer- 
ical and analytical approaches (including a mapping to a 
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FIG. 16: Total density versus time in a system of 500 lattice 
sites fully occupied at t = 0, reservoir densities g~ = 1, g + = 
0. Disordered case is average over 100 realizations drawn 
from a Gaussian distribution with fj = 2.0 and a v — 0.5. A 
disorder-induced slow-down in the decay to the steady state 
(g — 0.5) is clearly observed; presumably any oscillations in 
the approach to the steady state are averaged out. 



localization transition in an equivalent problem) to quan- 
tify these changes and shown that our results compare 
reasonably well with Monte Carlo simulations. 

While we believe that this mean-field discussion repro- 
duces qualitatively the effects of adding disorder, an ex- 
act treatment would be expected to provide better quan- 
titative agreement with simulations together with further 
physical understanding. Some progress has already been 
made (see the review by Stinchcombe [2^ ) . For example 
a form of the Harris Criterion [2j| can be applied to the 
ASEP which suggests that disorder should be relevant 
in the sense of introducing new critical behavior. And 
this new critical behavior can be elucidated by adding 



disorder to renormalization schemes developed for the 
pure case (see e.g., In carrying out such scaling 

of distributions one re-encounters many of the concepts 
highlighted above such as Griffith's phases and the im- 
portance of tails of the distribution. 

Many of the ideas discussed (e.g., Griffiths's phases, 
localization, importance of boundary conditions) might 
also be expected to apply to other nonequilibrium situa- 
tions and it would be interesting to see if we can apply 
generalizations of the methodology of Sections IIIIHVII to 
other problems. In particular we have studied a simple 
two-lane traffic model in which the fundamental dia- 
gram can have a double-peak structure (the ASEP with 
next-nearest neighbor interactions as studied by Popkov 
and Schiitz 2] also has such a double maxima). We would 
expect that disorder flattens the tops of these maxima 
leading to corresponding changes in the phase diagram 
but more interesting effects are also possible such as a 
relative change in the height of the two peaks. Studying 
quasi- ID models such as this two-lane system might also 
provide a bridge to understanding the effects of disorder 
on higher dimensional systems where one expects to find 
a wider range of possible disorder-induced effects. 

In conclusion, we hope that this mean-field treatment 
of the disordered asymmetric simple exclusion process 
provides a flavor of the general phenomena present in 
nonequilibrium models with quenched substitutional dis- 
order. There is much scope for further work on this and 
related models. 
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